The Role of Numerical Integration in Numerical Homogenization

Finite elements methods (FEMs) with numerical integration play a central role in numerical homogenization methods for partial differential equations with multiple scales, as the effective data in a homogenization problem can only be recovered from a microscopic solver at a finite number of points in the computational domain. In a multiscale framework the convergence of a FEM with numerical integration applied to the effective (homogenized) problem guarantees that the so-called macroscopic solver is consistent and convergent. Convergence results for FEM with numerical integration are however scarce in the literature and need often to be derived as a first step to analyze a numerical homogenization method for a given problem. In this paper we review and explain the main ideas in deriving convergence results for FEM with numerical integration for linear and nonlinear elliptic problems and explain the role of these methods in numerical homogenization. Résumé. Les méthodes d’éléments finis avec intégration numérique par quadrature jouent un rôle central dans l’homogénéisation numérique des équations aux dérivées partielles multi-échelles. En effet, les coefficients de l’équation homogénéisée ne peuvent être déterminés que pour un nombre fini de points du domaine considéré. Dans le cadre des méthodes multi-échelles basées sur un schéma macroscopique avec intégration numérique en la variable macroscopique et couplées à des schémas microscopiques autour des points de quadratures, la convergence d’une méthode d’élément fini avec intégration numérique pour le problème homogénéisé garantit que la méthode macroscopique est consistante et convergente. Comme les résultats de convergence pour les méthodes d’éléments finis avec intégration numérique ne sont connus que pour un nombre limité de problèmes, des nouveaux résultats de ce type sont alors un premier pas indispensable pour établir la convergence d’une méthode d’homogénéisation numérique pour un problème donné. Dans ce papier, nous effectuons un survol des idées clés pour l’analyse des méthodes d’éléments finis avec intégration numérique pour des problèmes linéaires et non linéaires et expliquons le rôle de ces méthodes dans l’homogénéisation numérique. Introduction When partial differential equations (PDEs) have coefficients that vary over multiple scales (highly oscillatory), classical numerical methods such as the finite element method (FEM), the finite volume method (FVM) or the finite difference method (FDM) are inefficient, as the mesh involved in any of these methods needs to resolve the smallest scale in the PDE to recover usual convergence rates. Such scale resolution is often prohibitive in many applications and there is a need for other type of numerical methods, called multiscale methods. ∗ This work is partially supported by the Swiss National Foundation under Grant 200021 134716/1. 1 ANMC, Section de Mathématiques, École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland, ([email protected]) c © EDP Sciences, SMAI 2014 2 ESAIM: PROCEEDINGS AND SURVEYS Among the large variety of multiscale problems, we concentrate on numerical methods for homogenization problems. Homogenization is a mathematical theory that study the convergence of PDEs with rapidly oscillating coefficients towards an averaged PDE describing the macroscopic or effective behavior of the physical process modeled by the differential equation [18,31,32]. Numerical homogenization methods are numerical methods able to approximate the effective solution of a highly oscillatory PDE without resolving the full oscillatory equation by direct discretization (see [3, 29] for recent reviews). Such methods are also called multiscale methods as they typically couple different solvers at different scales. The numerical analysis of numerical homogenization methods has been studied by many authors and we refer to [15] for a contribution of a pioneer in the field and to [3] for a recent review with numerous bibliographical entries. In this paper, we first discuss in Section 1 the typical structure of a numerical homogenization method. We then discuss the central role of FEM with numerical integration for such methods. We then review numerical homogenization methods for linear elliptic homogenization problems, nonlinear monotone elliptic homogenization problems and nonlinear nonmonotone elliptic homogenization problems in Sections 2,3,4, respectively. As our focus is on numerical integration we explain the various steps needed to control this error in a numerical homogenization method. We will show that for numerical homogenization, in order to obtain quantitative error bounds, FEM with numerical integration cannot be avoided as by its nature the macroscopic solver can only be defined at a finite number of points, which should be quadrature nodes of a macroscopic mesh. This raises the question of quantitative error bounds results for FEM with numerical integration for the various aforementioned problems. While the literature for the analysis of the FEM is abundant, works concerned with the analysis of the FEM with numerical integration are sparse even though for practical problems with non-constant parameters numerical integration is unavoidable. For linear elliptic problems, this problem has been thoroughly analyzed by Strang & Fix [37] and Ciarlet & Raviart [23]. For linear parabolic and hyperbolic problems results by Raviart [35] and Baker & Dougalis [17] are available. However for nonlinear problems quantitative error bounds have only recently be obtained for elliptic problems (of monotone and nonmonotone type) motivated by numerical homogenization [7, 11,12]. We close this introduction by mentioning that while we concentrate on a specific numerical homogenization method, namely the finite element heterogeneous multiscale method (FE-HMM) built in the framework of the heterogeneous multiscale method (HMM) [25], the link between numerical homogenization method and FEM with numerical integration is of general nature and not restricted to the specific framework described here. Notations. Let Ω ⊂ R be open and denote by W (Ω) the standard Sobolev spaces. We use the standard Sobolev norms ‖·‖W s,p(Ω) and semi-norms | · |W s,p(Ω). For p = 2 we will use the notations ‖·‖Hs(Ω) and | · |Hs(Ω), respectively. For k ∈ N and 1 ≤ p ≤ ∞, the broken norms ‖ · ‖W̄k,p(Ω) are given by ‖v‖W̄k,p(Ω) = (∑ K∈TH ‖v‖ p Wk,p(K) )1/p , if p <∞, ‖v‖W̄k,∞(Ω) = maxK∈TH ‖v‖Wk,∞(K), where TH is a partition of the closure of Ω, i.e., Ω = ⋃ K∈TH K, W (K) and ‖ · ‖Wk,p(K) denote the usual Sobolev spaces and norms on a closed subset K of Ω, respectively, and W̄ (Ω) is written as H̄(Ω). The notations |v|W̄k,p(Ω), |v|W̄k,∞(Ω) will be used for the corresponding semi-norm. 1. Numerical homogenization methods The general methodology of a numerical homogenization method such as the FE-HMM can be described as follows. Let Ω be an open bounded polygonal domain in R, V be a Hilbert space and consider the following (multiscale) problem: find u ∈ V such that L(u, a) = f in Ω, (1) ESAIM: PROCEEDINGS AND SURVEYS 3 with appropriate boundary conditions. Here L denotes a differential operator, a highly oscillatory data and f a right-hand side . The multiscale nature of the data, the operator and the solution is emphasized by the superscript ε (representing the typical size of a small scale in the considered problem). Assume now that u converges (weakly) in V to u as ε → 0, where the function u (a homogenized solution) solves a homogenized problem of the form L(u, a) = f in Ω. (2) In this latter problem, the small scales have been averaged out in the homogenization process. First, when one is not interested in the fine scale details involved in the problem (1) it is attractive to solve (2) numerically instead of (1), as (2) can in principle be solved with a standard FEM. Second, when fine scale features are required (often only in part of the computational domain) one can add corrector functions on the homogenized solution u to recover small scale information [34] or couple the solution of (1) in a part ω of the domain with the homogenized solution u in Ω\ω (in this direction see the recent results [9,16]). However, even when the primary goal is to solve numerically (2), this cannot be done straightforwardly as the effective data a(x), x ∈ Ω, are not known explicitly (except in very special cases) and must be computed with a micro FEM for a problem involving the original operator (1) on patches Kδj = xj + δY (sampling domains) around a given point xj ∈ Ω (here δ ≥ ε is a parameter usually of the size of ε and Y = (−1/2, 1/2)). Hence a(x) can only be approximated at a finite number of points xj , j = 1, . . . , N, in the computational domain. Let us describe a general framework for numerical homogenization methods such as the FE-HMM. Such methods rely on • a macroscopic domain Ω with a family of macroscopic triangulations TH such that ⋃ K∈TH K = Ω, where H is the maximum diameter of the element K ∈ TH and Ω is the closure of Ω; • a family of microscopic domains (sampling domains) Kδj = xj + δY, j = 1, . . . , N, and microscopic triangulation Th such that ⋃ T∈Th T = Kδj , where h is the maximum diameter of the element T ∈ Th. A numerical homogenization method (such as the FE-HMM) relies on (at least) two solvers, and on a data recovery process (1) a macroscopic solver LHMM for the effective problem L with a priori unknown data {a(xj)}j=1 defined on a macroscopic finite dimensional subspace VH(Ω) of V based on piecewise polynomials on each element K of the macroscopic triangulation TH of Ω; (2) a microscopic solver involving the operator L constrained by the macroscopic state (usually through boundary conditions involving the macro solution) defined on microscopic finite element spaces Vh(Kδj ) based on piecewise polynomials on each element T of the microscopic triangulation Thj ; (3) a data recovery process in which the effective data a(xj) at the point xj are computed using a suitable average involving the fine scale data a and the microscopic finite element solutions in each Kδj . We then consider the following problem: find uHMM ∈ VH such that LHMM (uHMM ) = f in Ω. (3) Fundamental questions are now the following. How should we choose the nodes xj in Ω to obtain the optimal accuracy of the numerical homogenization method with a minimal cost ? Recall that the computation of a(xj) for each xj is expensive as it involves a (micro) boundary value problem in a patch around xj . How should we choose the macro and micro meshsizes and the polynomial degree in order to obtain a given accuracy with the minimal computational cost ? We note that the macroscopic meshsize H is usually both larger and independent of ε and its size is solely dictated by the required accuracy for the effective problem. The microscopic meshsize is smaller than ε but only on patches Kδj that are usually also of size δ ' ε. A fully discrete a priori error estimates [1] then reveals that both meshes have to be refined simultaneously and that the computational cost, independent of ε, is proportional to the macroscopic degrees of freedom. 1For simplicity we do not consider a highly oscillatory right-hand side (e.g., fε) but we note that this situation can usually be handled with minor modifications of the method described below (see e.g., [10]). 4 ESAIM: PROCEEDINGS AND SURVEYS For the analysis of such a method, we usually introduce an intermediate problem, namely a macroscopic solver LH for the effective problem L with data {a(xj)}j=1 assumed to be exact and computed at the same location xj ∈ Ω as for the FE-HMM. We thus have three macroscopic solutions: u the solution of (2), uHMM the solution of (3) and uH the solution of LH(uH) = f in Ω. The analysis of the FE-HMM relies on the following decomposition [4] ‖u− uHMM‖ ≤ ‖u− uH‖+ ‖uH − uHMM‖. (4) The first term in the right-hand side of (4) is usually called the macroscopic error. It depends on the type of discretization used at the macroscopic level, the macroscopic meshsize and on the quadrature points xj chosen in the computational domain to recover {a(xj)}. The second term in the right-hand side of (4) is usually called the macroscopic error comprise the so-called microscopic and modeling errors and depends on the type of discretization used at the microscopic level, the macroscopic meshsize, the sampling domain size and the type of boundary conditions used for the micro solver. We note that it is the second term that involves the mathematical homogenization theory for its analysis, while the control of the first term involves a classical question in the numerical analysis of the FEM as described in the introduction, which goes back to Strang & Fix [37] and Ciarlet & Raviart [23] for linear elliptic problem. For numerical homogenization methods, as we will show in this paper, FEM with numerical integration cannot be avoided as by its nature the macroscopic solver can only be defined at quadrature nodes. Thus any numerical homogenization method for a new class of problems must come with its companion analysis of FEM with numerical integration. These results are however not trivial already in the linear case (see [23] for a very general analysis). In particular for nonlinear problems, quantitative error bounds have only recently be obtained for elliptic problems motivated by numerical homogenization [7, 11, 12]. We will highlight in the following sections the main ideas to derive error estimates for FEM with numerical integration for various linear and nonlinear problems and show how such solvers naturally arise in numerical homogenization. We close this introduction by noting that while we only describe here multiscale methods coupling elliptic equations at both the macro and the micro scales, one could also couple different physics at different scales. This has been pursued for example in [6], where the macroscopic solver discretizes an effective Darcy equation, while the microscopic solver discretizes Stokes problems around quadrature points of the macromesh. There again, results for FEM with numerical integration for the Darcy problem are important. 2. Numerical homogenization and numerical integration for linear elliptic problems In this section we first detail the abstract method described in the introduction for elliptic linear multiscale problems of the form 2 L(u, a) = −∇ · (a∇u) = f in Ω, u = 0 on ∂Ω. (5) We assume that f ∈ L(Ω) and that the family of tensors a(x) ∈ (L∞(Ω))d×d indexed by ε > 0 is uniformly elliptic and bounded, thus for any ε there exist a unique solution u ∈ H 0 (Ω) of Problem (5) and the family of solution {u} is bounded independently of ε. For such problems the homogenization theory [18, 31, 32] ensures the existence of a subsequence of {u} that weakly converges in H 0 (Ω) to function u, solution of the homogenized problem L(u, a) = −∇ · (a∇u) = f in Ω with u = 0 on ∂Ω that reads in weak form: find u ∈ H 0 (Ω) such that